R Notebook

Markdown instruction

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

Import Library

#install.packages("dplyr")
library(ggplot2)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(corrplot)

corrplot 0.92 loaded

library(readxl)
library(lmtest) 

Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

library(forecast)
library(DIMORA)

Loading required package: minpack.lm
Loading required package: numDeriv
Loading required package: reshape2
Loading required package: deSolve

library(fpp2)

── Attaching packages ────────────────────────────────────────────────────────────────────────── fpp2 2.5 ──
✔ fma       2.5     ✔ expsmooth 2.3

Data upload

btc_price <- read.csv("BTC_price_Bitcoin.csv", sep=",", header = TRUE) # daily
dollar_circulation <- read.csv("CURRCIR.csv", sep=",", header = TRUE) # monthly
btc_addr_total <- read.csv("BTC_Total Addresses_Bitcoin.csv", sep=",", header = TRUE) # daily
btc_hashrate_difficulty <- read.csv("BTC_Hash Rate_Difficulty.csv", sep=",", header = TRUE) # daily
btc_transaction <- read.csv("BTC_Number of Transactions_Bitcoin.csv", sep=",", header = TRUE) # daily
btc_miner_rewards <- read.csv("BTC_Miner Rewards_undefined.csv", sep=",", header = TRUE) # daily

btc_avg_fees <- read.csv("BTC_Average Transaction Fees_Bitcoin.csv", sep=",", header = TRUE) # daily not null from 29/11/09
btc_days_halving <- read.csv("days_to_halving.csv", sep=",", header = TRUE) # daily

Function

# Normalize data
norm <- function(lista_numeri) {
  # Calcola il minimo e il massimo della lista
  minimo <- min(lista_numeri)
  massimo <- max(lista_numeri)
  
  # Normalizza la lista tra 0 e 1
  lista_normalizzata <- (lista_numeri - minimo) / (massimo - minimo)
  
  # Restituisci la lista normalizzata
  return(lista_normalizzata)
}
# From daily to monthly
monthly <- function(data, col_date) {
  # Assicurati che la colonna 'DateTime' sia di tipo Date
  data[[col_date]] <- as.Date(data[[col_date]], format = "%Y-%m-%dT%H:%M:%S.000Z")
  
  # Estrai il mese dalla colonna 'DateTime'
  data <- data %>%
    mutate(Month = format(get(col_date), "%Y-%m"))
  
  # Seleziona solo colonne numeriche per calcolare la media
  cols_num <- names(data)[sapply(data, is.numeric)]
  
  # Calcola la media delle colonne numeriche per ogni mese
  new_dataset <- data %>%
    group_by(Month) %>%
    summarise(across(all_of(cols_num), mean))
  
  new_dataset$Month <- as.Date(paste(new_dataset$Month, "01", sep = "-"))
  
  # Restituisci il nuovo dataset
  return(new_dataset)
}

EDA

From day to month

btc_price <- monthly(btc_price, "DateTime")
btc_addr_total <- monthly(btc_addr_total, "DateTime")
btc_hashrate_difficulty <- monthly(btc_hashrate_difficulty, "DateTime")
btc_transaction <- monthly(btc_transaction, "DateTime")
btc_miner_rewards <- monthly(btc_miner_rewards, "DateTime")
btc_avg_fees <- monthly(btc_avg_fees, "DateTime")
#btc_days_halving <- monthly(btc_days_halving, "date")

Global date

start_date <- as.Date("2011-09-01") # 2011-09-01
end_date <- as.Date("2023-10-01") # 2023-10-01

Date

btc_price <- btc_price[btc_price$Month >= start_date & btc_price$Month <= end_date, ] # price
dollar_circulation <- dollar_circulation[dollar_circulation$DATE >= start_date & dollar_circulation$DATE <= end_date, ] # price
btc_addr_total <- btc_addr_total[btc_addr_total$Month >= start_date & btc_addr_total$Month <= end_date, ] # price
btc_hashrate_difficulty <- btc_hashrate_difficulty[btc_hashrate_difficulty$Month >= start_date & btc_hashrate_difficulty$Month <= end_date, ] # price
btc_transaction <- btc_transaction[btc_transaction$Month >= start_date & btc_transaction$Month <= end_date, ] # price

btc_miner_rewards <- btc_miner_rewards[btc_miner_rewards$Month >= start_date & btc_miner_rewards$Month <= end_date, ] # price

Error in btc_miner_rewards

nrows <- nrow(btc_miner_rewards)
row_pre <- btc_miner_rewards[(nrows - 29):(nrows - 15), -1]
btc_miner_rewards[(nrows - 14):nrows, -1] <- row_pre

Log Data

btc_price$log_Price <- log(btc_price$Price)
dollar_circulation$log_CURRCIR <- log(dollar_circulation$CURRCIR)
btc_addr_total$log_total <- log(btc_addr_total$Total.With.Balance)
btc_hashrate_difficulty$log_HashRate <- log(btc_hashrate_difficulty$Hash.Rate)
btc_transaction$log_transaction <- log(btc_transaction$Number.Of.Transactions)
btc_miner_rewards$log_Rewards <- log(btc_miner_rewards$Rewards)

Plot

ggplot(btc_price, aes(x = Month, y = log_Price)) +
  geom_line(color = "orange", size = 1) +
  labs(x = "Date", y = "Price Bitcoin", title = "Bitcoin price over time")

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

Plot

ggplot() +
  geom_line(data = btc_price, aes(x = Month, y = norm(log_Price)), color = "orange", size = 1, linetype = "solid") +
  geom_line(data = dollar_circulation, aes(x = as.Date(DATE), y = norm(log_CURRCIR)), color = "red", size = 0.2, linetype = "solid") +
  geom_line(data = btc_addr_total, aes(x = Month, y = norm(log_total)), color = "green", size = 0.2, linetype = "solid") +
  geom_line(data = btc_transaction, aes(x = Month, y = norm(log_transaction)), color = "black", size = 0.2, linetype = "solid") +
  geom_line(data = btc_hashrate_difficulty, aes(x = Month, y = norm(log_HashRate)), color = "blue", size = 0.2, linetype = "solid") +
  geom_line(data = btc_miner_rewards, aes(x = Month, y = norm(log_Rewards)), color = "blue", size = 0.2, linetype = "solid") +

  labs(x = "Date", y = "Log(Price)", title = "Bitcoin Price and other Time Series") +
  scale_linetype_manual(name = "Legend", values = c("solid", "dashed"), labels = c("Bitcoin Price", "Another Time Series"))

Create Dataset

dataset_log <- data.frame(price = btc_price$log_Price, 
                      dollar_circulation = dollar_circulation$log_CURRCIR,
                      total_addr = btc_addr_total$log_total,
                      transaction = btc_transaction$log_transaction,
                      hashrate = btc_hashrate_difficulty$log_HashRate,
                      rewards = btc_miner_rewards$log_Rewards)

dataset <- data.frame(price = btc_price$Price, 
                      dollar_circulation = dollar_circulation$CURRCIR,
                      total_addr = btc_addr_total$Total.With.Balance,
                      transaction = btc_transaction$Number.Of.Transactions,
                      hashrate = btc_hashrate_difficulty$Hash.Rate,
                      rewards = btc_miner_rewards$Rewards)

acf(dataset$price)

mat <- cor(dataset)

corrplot(mat, method = "color", addCoef.col = "black")

Differentation

dataset$price_diff <- c(NA, diff(dataset$price))
dataset$dollar_circulation_diff <- c(NA, diff(dataset$dollar_circulation))
dataset$total_addr_diff <- c(NA, diff(dataset$total_addr))
dataset$transaction_diff <- c(NA, diff(dataset$transaction))
dataset$hashrate_diff <- c(NA, diff(dataset$hashrate))
dataset$rewards_diff <- c(NA, diff(dataset$rewards))

# Rimuovi le variabili temporali originali
dataset_senza_tempo <- subset(dataset, select = -c(price, dollar_circulation, total_addr, transaction, hashrate, rewards))
# Calcola la matrice di correlazione
matrix_correlazione <- cor(dataset_senza_tempo, use = "complete.obs")
# Visualizza la matrice di correlazione
print(matrix_correlazione)

                        price_diff dollar_circulation_diff total_addr_diff transaction_diff hashrate_diff
price_diff              1.00000000              0.06055474     0.229407118       0.15074382    0.29163970
dollar_circulation_diff 0.06055474              1.00000000     0.101964436      -0.04758256   -0.11206328
total_addr_diff         0.22940712              0.10196444     1.000000000       0.20045376    0.06110231
transaction_diff        0.15074382             -0.04758256     0.200453760       1.00000000    0.03003450
hashrate_diff           0.29163970             -0.11206328     0.061102310       0.03003450    1.00000000
rewards_diff            0.01177412             -0.06658894    -0.002962721       0.01710025    0.03394439
                        rewards_diff
price_diff               0.011774119
dollar_circulation_diff -0.066588944
total_addr_diff         -0.002962721
transaction_diff         0.017100247
hashrate_diff            0.033944391
rewards_diff             1.000000000

corrplot(matrix_correlazione, method = "color", addCoef.col = "black")

Linear Regression

price <- log(dataset$price)
tt <- 1:NROW(dataset)

plot(tt, price, xlab="Time", ylab="Bitcoin price")

autocorrelation function

acf(price)

fit a linear regression model

fit1 <- lm(price~ tt)
summary(fit1)

Call:
lm(formula = price ~ tt)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.11319 -0.35944  0.01361  0.64569  2.17326 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.923081   0.159401   18.34   <2e-16 ***
tt          0.060502   0.001881   32.16   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9581 on 144 degrees of freedom
Multiple R-squared:  0.8778,    Adjusted R-squared:  0.8769 
F-statistic:  1034 on 1 and 144 DF,  p-value: < 2.2e-16

plot(tt, price, xlab="Time", ylab="Bitcoin price")
abline(fit1, col=3)

check the residuals? are they autocorrelated? Test of DW

dwtest(fit1)

    Durbin-Watson test

data:  fit1
DW = 0.062633, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0

check the residuals

resfit1<- residuals(fit1)
plot(resfit1,xlab="Time", ylab="residuals" )

let us do the same with a linear model for time series, so we transform the data into a ‘ts’ object

price.ts <- ts(price, frequency = 4)
ts.plot(price.ts, type="o")

## we fit a linear model with the tslm function
fitts<- tslm(price.ts~trend)

###obviously it gives the same results of the first model
summary(fitts)

Call:
tslm(formula = price.ts ~ trend)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.11319 -0.35944  0.01361  0.64569  2.17326 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.923081   0.159401   18.34   <2e-16 ***
trend       0.060502   0.001881   32.16   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9581 on 144 degrees of freedom
Multiple R-squared:  0.8778,    Adjusted R-squared:  0.8769 
F-statistic:  1034 on 1 and 144 DF,  p-value: < 2.2e-16

dwtest(fitts)

    Durbin-Watson test

data:  fitts
DW = 0.062633, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0

Linear regression with trend and seasonality and forecasting exercise

#take a portion of data and fit a linear model with tslm
price1 <- window(price.ts, start=1, end=31 -.1)
plot(price1)

m1<- tslm(price1 ~ trend+season)
summary(m1)

Call:
tslm(formula = price1 ~ trend + season)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6729 -0.5493 -0.1696  0.5580  2.3361 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.464363   0.209349  11.772   <2e-16 ***
trend        0.070706   0.002301  30.725   <2e-16 ***
season2     -0.043150   0.225364  -0.191    0.848    
season3     -0.012176   0.225400  -0.054    0.957    
season4      0.033072   0.225458   0.147    0.884    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8728 on 115 degrees of freedom
Multiple R-squared:  0.8916,    Adjusted R-squared:  0.8878 
F-statistic: 236.4 on 4 and 115 DF,  p-value: < 2.2e-16

fit<- fitted(m1)

plot(price1)
lines(fitted(m1), col=2)

fore <- forecast(m1)
plot(fore)

#analysis of residuals
res<- residuals(m1) 
plot(res) 

#the form of residuals seems to indicate the presence of negative autocorrelation
Acf(res)

dw<- dwtest(m1, alt="two.sided")
dw

    Durbin-Watson test

data:  m1
DW = 0.086005, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is not 0

Nonlinear models for new product growth (diffusion models)

bm_tw<-BM(price,display = T)

summary(bm_tw)

Call: ( Standard Bass Model )

  BM(series = price, display = T)

Residuals:
     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-21.2146  -5.5677   0.4326  -1.7988   3.5383   7.2830 

Coefficients:
 

---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Residual standard error  7.767938  on  143  degrees of freedom
 Multiple R-squared:   0.999796  Residual sum of squares:  8628.743

pred_bmtw<- predict(bm_tw, newx=c(1:146))
pred.insttw<- make.instantaneous(pred_bmtw)


plot(price, type= "b",xlab="Quarter", ylab="Quarterly revenues",  pch=16, lty=3, cex=0.6, xlim=c(1,146))
lines(pred.insttw, lwd=2, col=2)

###GBMr1
GBMr1tw<- GBM(price,shock = "rett",nshock = 1,prelimestimates = c(4.463368e+04, 1.923560e-03, 9.142022e-02, 24,38,-0.1))

######GBMe1

GBMe1tw<- GBM(price,shock = "exp",nshock = 1,prelimestimates = c(4.463368e+04, 1.923560e-03, 9.142022e-02, 12,0.1,0.1))

summary(GBMe1tw)

Call: ( Generalized Bass model with 1  Exponential  shock )

  GBM(series = price, shock = "exp", nshock = 1, prelimestimates = c(44633.68, 
    0.00192356, 0.09142022, 12, 0.1, 0.1))

Residuals:
     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-5.83499 -2.24663 -0.08100 -0.08908  1.91749  6.93016 

Coefficients:
 

---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Residual standard error  3.165267  on  140  degrees of freedom
 Multiple R-squared:   0.99991  Residual sum of squares:  1402.648

pred_GBMe1tw<- predict(GBMe1tw, newx=c(1:146))
pred_GBMe1tw.inst<- make.instantaneous(pred_GBMe1tw)

plot(price, type= "b",xlab="Quarter", ylab="Quarterly revenues",  pch=16, lty=3, cex=0.6, xlim=c(1,146))
lines(pred_GBMe1tw.inst, lwd=2, col=2)

######GGM 
GGM_tw<- GGM(price, prelimestimates=c(4.463368e+04, 0.001, 0.01, 1.923560e-03, 9.142022e-02))

Warning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs produced

summary(GGM_tw)

Call: ( Guseo Guidolin Model )

  GGM(series = price, prelimestimates = c(44633.68, 0.001, 0.01, 
    0.00192356, 0.09142022))

Residuals:
    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-7.3569 -2.6948  0.5882  0.2721  3.2064  8.1475 

Coefficients:
 

---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Residual standard error  4.118758  on  141  degrees of freedom
 Multiple R-squared:   0.999847  Residual sum of squares:  2391.947

pred_GGM_tw<- predict(GGM_tw, newx=c(1:146))
pred_GGM_tw.inst<- make.instantaneous(pred_GGM_tw)

plot(price, type= "b",xlab="Quarter", ylab="Quarterly revenues",  pch=16, lty=3, cex=0.6, xlim=c(1,146))
lines(pred_GGM_tw.inst, lwd=2, col=2)

lines(pred.insttw, lwd=2, col=3)

###Analysis of residuals
res_GGMtw<- residuals(GGM_tw)
acf<- acf(residuals(GGM_tw))

fit_GGMtw<- fitted(GGM_tw)
fit_GGMtw_inst<- make.instantaneous(fit_GGMtw)

####SARMAX refining
library(forecast)
####SARMAX model with external covariate 'fit_GGM' 
s2 <- Arima(cumsum(price), order = c(3,0,1), seasonal=list(order=c(3,0,1), period=4),xreg = fit_GGMtw)
summary(s2)

Series: cumsum(price) 
Regression with ARIMA(3,0,1)(3,0,1)[4] errors 

Coefficients:
         ar1      ar2     ar3     ma1     sar1    sar2    sar3    sma1  intercept    xreg
      1.9142  -0.9560  0.0167  0.1962  -0.8595  0.0826  0.0883  0.7991    -0.9844  1.0030
s.e.  0.0203   0.0407  0.0220  0.0558   0.0118  0.0050  0.0057  0.0392     1.3369  0.0024

sigma^2 = 0.04182:  log likelihood = 25.57
AIC=-29.14   AICc=-27.17   BIC=3.67

Training set error measures:
                       ME      RMSE       MAE         MPE      MAPE      MASE        ACF1
Training set 0.0009240674 0.1973691 0.1483844 -0.01038485 0.3713901 0.0200291 0.001633754

pres2 <- make.instantaneous(fitted(s2))


plot(price, type= "b",xlab="Quarter", ylab="Quarterly revenues",  pch=16, lty=3, xaxt="n", cex=0.6)
lines(fit_GGMtw_inst, lwd=1, lty=2)

lines(pres2, lty=1,lwd=1)

ARIMA models

library(fpp2)

── Attaching packages ────────────────────────────────────────────────────────────────────────── fpp2 2.5 ──
✔ ggplot2   3.4.2      ✔ fma       2.5   
✔ forecast  8.21.1     ✔ expsmooth 2.3   

library(forecast) 
?fpp2

plot(price)

Acf(price)

Pacf(price)

pricets<- tsdisplay(price)

###General indication: if the ACF is exponentially decaying or sinusoidal and there is a significant spike at lag p in PACF and nothing else, 
##it may be an ARMA(p,d,0). If the PACF is exponentially decaying or sinusoidal and there is a significant spike at lag p in ACF and nothing else, it may be an ARMA(0,d,q). 

arima1<- Arima(price, order=c(0,0,3))
summary(arima1)

Series: price 
ARIMA(0,0,3) with non-zero mean 

Coefficients:
         ma1     ma2     ma3    mean
      1.7886  1.6759  0.8049  7.3514
s.e.  0.0556  0.0583  0.0531  0.2606

sigma^2 = 0.3742:  log likelihood = -135.94
AIC=281.89   AICc=282.32   BIC=296.81

Training set error measures:
                    ME      RMSE       MAE       MPE     MAPE     MASE      ACF1
Training set 0.0147935 0.6032425 0.4836402 -6.579356 11.83727 2.842824 0.5551343

resid1<- residuals(arima1)
tsdisplay(resid1)

plot(price)
lines(fitted(arima1), col=2)

for1<- forecast(arima1)
plot(for1)

First Arima model

plot(price, ylab="retail index",xlab="year")

tsdisplay(price)

##first difference
diff1<- diff(price) 
###seasonal difference
diff4<- diff(price, lag=4) 
tsdisplay(diff1)

tsdisplay(diff4)

####first Arima model 
a1<- Arima(price, order=c(0,1,1), seasonal=c(0,0,1))
fit1<- fitted(a1)

plot(price)
lines(fit1, col=2)

f1<- forecast(a1)
plot(f1)

r1<- residuals(a1)
tsdisplay(r1) 

Second Arima model

a2<- Arima(price, order=c(0,1,1), seasonal=c(0,0,2))
fit2<- fitted(a2)

plot(price)
lines(fit2, col=2)

f2<- forecast(a2)
plot(f2)

r2<- residuals(a2)
tsdisplay(r2) 

Third Arima model

a3<- Arima(price, order=c(0,1,1), seasonal=c(0,1,1))
fit3<- fitted(a3)

plot(price)
lines(fit3, col=2)

f3<- forecast(a3)
plot(f3)

r3<- residuals(a3)
tsdisplay(r3) 

Fourth Arima model

a4<- Arima(price, order=c(0,1,2), seasonal=c(0,1,1))
fit4<- fitted(a4)

plot(price)
lines(fit4, col=2)

f4<- forecast(a4)
autoplot(f4)

r4<- residuals(a4)
tsdisplay(r4) 

Fifth Arima model

auto.a<- auto.arima(price)
auto.a

Series: price 
ARIMA(1,1,0) with drift 

Coefficients:
         ar1   drift
      0.3878  0.0566
s.e.  0.0778  0.0298

sigma^2 = 0.04939:  log likelihood = 13.26
AIC=-20.52   AICc=-20.35   BIC=-11.59

autoplot(forecast(auto.a))

checkresiduals(auto.a)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,0) with drift
Q* = 22.133, df = 9, p-value = 0.008466

Model df: 1.   Total lags used: 10

##plot of seasonal differentiated data, with ACF and PACF
tsdisplay(diff(price,12), main="seasonal differenced data", xlab="year")

##We fit an ARIMA model base on the inspection of ACF and PACF with 3 AR components
fit<- Arima(price, order=c(3,0,1),seasonal=c(0,1,2), lambda=0)
summary(fit)

Series: price 
ARIMA(3,0,1) with non-zero mean 
Box Cox transformation: lambda= 0 

Coefficients:
         ar1      ar2     ar3     ma1    mean
      1.4137  -0.7747  0.3576  0.1445  1.5330
s.e.  0.1839   0.2713  0.1148  0.1833  0.7488

sigma^2 = 0.003475:  log likelihood = 205.78
AIC=-399.55   AICc=-398.95   BIC=-381.65

Training set error measures:
                     ME      RMSE       MAE       MPE     MAPE    MASE       ACF1
Training set 0.05831207 0.2394853 0.1746367 0.8848972 3.383428 1.02651 -0.2082698

##Check residuals
tsdisplay(residuals(fit))

Box.test(residuals(fit), lag=36,fitdf=6, type="Ljung")

    Box-Ljung test

data:  residuals(fit)
X-squared = 37.175, df = 30, p-value = 0.1721

##perform the forecasting
f<- forecast(fit)

plot(f, ylab="sales", xlab="year")

air.model <- Arima(window(price,end=146),order=c(0,1,1),
                   seasonal=list(order=c(0,1,1),period=12),lambda=0)
plot(forecast(air.model,h=24))
lines(price)

# Apply fitted model to later data
air.model2 <- Arima(window(price,start=1),model=air.model)

# in-sample one-step forecasts
accuracy(air.model)

                      ME      RMSE       MAE        MPE     MAPE     MASE        ACF1
Training set -0.05360399 0.3327414 0.2287816 -0.6502599 3.658534 1.344772 -0.07160046

# out-of-sample one-step forecasts
accuracy(air.model2)

                      ME      RMSE       MAE        MPE     MAPE     MASE        ACF1
Training set -0.05360399 0.3327414 0.2287816 -0.6502599 3.658534 1.344772 -0.07160046

ARMAX Models

armax1<- Arima(price, xreg=dataset_log$dollar_circulation, order=c(1,0,1))
res1<- residuals(armax1)
Acf(res1)

fitted(armax1)

Time Series:
Start = 1 
End = 146 
Frequency = 1 
  [1]  2.091489  1.825582  1.317792  1.079898  1.426339  2.108251  1.679215  1.690872  1.700841  1.731494
 [11]  1.948983  2.201273  2.613059  2.521132  2.626462  2.507943  2.716105  2.849346  3.511862  4.258859
 [21]  5.105274  4.706627  4.698566  4.473554  4.910694  4.952853  5.294621  6.812591  6.613939  6.882095
 [31]  6.082504  6.412056  6.041619  6.232432  6.472638  6.353501  6.222363  6.005696  5.925920  6.000025
 [41]  5.787144  5.436091  5.561916  5.628957  5.423743  5.506027  5.509012  5.728103  5.498019  5.514887
 [51]  5.716613  6.000273  6.108872  6.002581  6.110043  6.053700  6.129470  6.217940  6.607589  6.441679
 [61]  6.372712  6.456957  6.567509  6.666671  6.771188  6.829174  7.180119  7.005705  7.173147  7.720319
 [71]  7.926871  7.805486  8.523932  8.191952  8.763107  9.081779  9.755043  9.236725  9.165211  8.989041
 [81]  9.023974  8.983845  8.728155  8.881322  8.787780  8.763677  8.792990  8.461557  8.134137  8.192134
 [91]  8.275328  8.302494  8.652074  9.025397  9.186619  9.253105  9.284655  9.101616  9.075293  9.000331
[101]  8.845539  9.082660  9.314844  8.943486  9.103742  9.258128  9.212633  9.298374  9.462879  9.266633
[111]  9.518255  9.878156 10.203930 10.539189 10.822939 11.065714 10.913371 10.636783 10.387967 10.451131
[121] 10.825671 10.704060 11.053093 11.011211 10.749139 10.555663 10.722945 10.652860 10.627915 10.297995
[131] 10.056031  9.952615 10.083163  9.867203  9.948428  9.772788  9.772431  9.979983 10.157293 10.189559
[141] 10.353882 10.233927 10.242435 10.315912 10.217728 10.169996

plot(price)
lines(fitted(armax1), col=2)